function [CT,PD,PT,OTHERS] = dgcomp_postp(CD,base,ioref,ref,phc)
% [CT,PD,PT,OTHERS] = dgcomp_postp(CD,base,ioref,ref,phc)
%
% Post-processing of the compressible DG solution: given the
% deviations from reference for the conservative variables CD,
% compute:
%  CT: conservative (rho,E,U), total
%  PD: primitive (p,T,u), deviations
%  PT: primitive (p,T,u), total
%  OTHERS: other useful variables (total): (pi,theta)
%
% Input values:
%  CD:    conservative, deviations
%  base:  finite element basis
%  ioref: 4D array with the output reference state
%  ref:   structure with the reference atmosphere
%  phc:   physical constants
%
% Notice that all the variables are expressed with respect to the
% basis functions.
%
% Notice that, is the reference atmosphere is not well defined, PD
% will be NaN. The remaining variables are still meaningful.

 % 1) Compute CT
 Cref = ioref;
 Cref(1,:,:) = Cref(1,:,:) + permute(ref.rho,[3 1 2]);
 Cref(2,:,:) = Cref(2,:,:) + permute(ref.e,  [3 1 2]);
 CT = CD + Cref;
 
 % 2) Compute PT and OTHERS
 [PT,OTHERS] = C2P(CT,base,ref,phc);
 
 % 3) Compute PD
 PD = PT - C2P(Cref,base,ref,phc);

return


function [P,PITHETA] = C2P(C,base,ref,phc)
% Conservative to primitive conversion: due to the nonlinear relations
% between conservative and primitive variables, the conversion
% requires two steps: variable transformation and projection on the
% basis function.

 % precompute the projection matrix
 % scaled basis functions
 for i=1:base.pk
   wgp(i,:) = base.wg .* base.p(i,:);
 end
 % mass matrix
 for i=1:base.pk
   for j=1:base.pk
     M(i,j) = sum( wgp(i,:) .* base.p(j,:) );
   end
 end
 MI = M^(-1);
 PL2 = MI * wgp;

 P = zeros(size(C));
 PITHETA = zeros(2,size(P,2),size(P,3));
 for ie=1:size(C,3) % loop on the elements

   % precompute the gravitational potential
   phi = ref.phi(:,ie)' * base.p;

   % conservative variables at the quadrature nodes
   rho = C(1,:,ie) * base.p;
   E   = C(2,:,ie) * base.p;
   for i=3:size(C,1)
     U(i,:) = C(i,:,ie) * base.p;
   end

   % primitive variables at the quadrature nodes
   % velocity
   for i=3:size(C,1)
     u(i,:) = U(i,:)./rho;
   end
   % temperature
   T = 1/phc.cv*( E./rho - 0.5*sum(u.^2,1) - phi );
   % pressure
   p = rho .* phc.rgas .* T;

   % pi
   pi = (p/phc.p_s).^phc.kappa;
   % theta
   theta = T./pi;

   % L2 projection
   P(1,:,ie) = PL2 * p';
   P(2,:,ie) = PL2 * T';
   for i=3:size(C,1)
     P(i,:,ie) = PL2 * u(i,:)';
   end

   PITHETA(1,:,ie) = PL2 * pi';
   PITHETA(2,:,ie) = PL2 * theta';

 end

return

